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Introduction 


The Langley Aerotherinodyuainic Upwind Relaxation Algorithm (LAURA) was originally de- 
veloped by Gnoiro 1 to solve steady-llow problems. The desire to validate the algorithm with shock 
tube experimental data 2 has motivated the development of a lime-accurate version of the LAURA 
code. Using the LAURA scheme in this manner is not the optimal approach to the solution of 
unsteady flow problems, but the code is easily modified to produce time-accurate results. 

The objective of this work is to develop and test the time-accurate version of LAURA. The 
algorithm’s modifications are discussed in the next section. Testing the solution accuracy of LAURA 
is accomplished by comparing the computational results with an exact solution. For this purpose 
the present study is solely concerned with inviscid, perfect-gas shock tube flow. Real gas effects 
will be examined in subsequent studies. The approach to analyzing the overall solution accuracy 
of LAURA is discussed in the section entitled Parameters of the Study. 


Time- Accurate Relaxation Algorithm 


The governing equation for inviscid flow of a perfect gas may be written 
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In Eq. 1 the first term describes the time rate of change of conserved quantity q in the control 
volume Q and the second term describes the convective flux f through the cell walls; n is the unit 
vector normal to the cell walls and a is the cell wall area. The finite-volume formulation of Eq. 1 
is given by 

["77"k + E [^+1 1 n/+i<7/+i - ft • n t <r, ] = 0 (2) 

Note in Eq. 2 that uppercase integer variables L denote computational coordinates at cell centers, 
and lowercase integer variables / denote cell faces. We define g/ = f/ • ni so that Eq. 2 becomes 

[%f-h + ZJ = o ( 3 ) 

/=t J,k 

where 

P 

pu 

q = po (4) 

pw 

. pE . 
pU 

pU u + pn x 

g = pUv + pn y (5) 

pUw + pn z 
pUJI 

In Eqs. 4 and 5 p represents density, u, v, and w are the scalar components of velocity, E is the 
total energy per unit mass, U is the normal component of velocity through the cell wall, P is the 
pressure, and II represents the total enthalpy per unit mass. 

The governing relaxation equation in the LAURA 1 code for inviscid, steady flow is defined by 

q l +1 = qZ + M^rj, (6) 

where Mj, is the point-implicit Jacobian given by 

Ml = ^ E [|A. + .h + . + |AiM" (?) 

and r is the right-hand-side solution (residual) vector given by 

r L = ~ ^2 [gJ+i^z+i “ g i c t\ n ( 8 ) 
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Spatial differencing of the residual vector can be either first or second order accurate. In Eq. 7 
rfiNV represents the inviscid relaxation factor. The matrix |A,| is the Roe’s averaged 3 Jacobian 
of g with respect to q at cell wall / with absolute, limited eigenvalues as described in Ref. 1. 

We employ M local relaxation steps before advancing a time step At to obtain first order 

accuracy ill time in the following manner. 

q n,m+i _ q n,m + Ml l r L (9) 

Mt = r -^p- £ t|A, + ,k, + l + IAiN ”’ 1 + <“» 

1 

and 

n "' m _ n "’ 1 

n = - £ [»+.»,+, - wr- 1 Hi) 

In Eq. 10 I represents the identity matrix. After M local iterations, the solution is advanced to 
time step n + 1 by setting 


<&*'•' = <if < i2 > 

The number of local iterations M should be large enough such that the residual error defined 
by Eq. 11 above is smaller than O( Al'). Extensions of this method to viscous, nonequilibrium flows 
are trivial, but as yet untested. 

Parameters of the Study 

The parametric study is conducted to define the uncertainty and computational cost associated 
with applying LAURA to unsteady flow problems. The parameters examined include Courant 
number, grid resolution, relaxation sweeps, and the inviscid relaxation factor. 

For the simple shock tube, Courant number defines the number of cells the shock wave prop- 
agates every time step. 

cn = ^ (W) 

OX 

where, A is the velocity of the shock wave, 6t is the time step, and 6x is the grid spacing. If Sx 
is constant then Courant number varies with St. The effect of the time step on a solution may 
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be examined by varying the Courant number. Decreasing the time step improves accuracy at the 

price of additional computational cost. 

The grid spacing is related to the number of cells in the grid. 

6x = Srid extent ( 14 ) 

cells 

Increasing the number of grid cells improves the grid resolution and solution accuracy. This im- 
proved accuracy carries with it additional computational cost. 

The number of relaxation sweeps is the number of iterations performed to “relax’ the solution 
ahead one time step. With a sufficient number of relaxation sweeps the solution converges to 
time-accuracy. There is a direct trade-off between relaxation sweeps and computational cost. This 
trade-off prompts an investigation into the minimum number of relaxation sweeps required to obtain 
a time-accurate converged solution. 

The in viscid relaxation factor is the only parameter to be studied with no associated trade-off in 
computational cost. The rfi NV was initially implemented in the algorithm to increase the stability 
of the scheme. Two values of rf tN v are examined in this study: rfi N v= 1-5 and rf INV = 2.0. 
Indirectly, this factor controls the speed with which a solution converges. 

The parametric study is conducted to analyze the capabilities of the time-accurate LAURA 
code in predicting unsteady How solutions. The four coupled parameters are isolated to examine 
the contribution of each to the overall accuracy of a LAURA solution. The study examines discrete 
values for each of the parameters. The results suggest preferred values for each parameter to 
efficiently achieve time-accurate solutions. These suggested values for the parameters, while not 
completely optimized, are selected by examining trends in the parametric study. 

Test Case and Computational Mesh 

The initial conditions of the simple shock tube test case are taken from Ref. 4. The low 
pressure side of the diaphragm is initialized with standard atmospheric conditions at sea level. The 
pressure ratio across the diaphragm is 10. The density ratio is 8. Both sides have zero initial 
velocity. The shock tube test case assumes inviscid, perfect gas (7 = 1.4). 
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A one-dimensional computational grid extending from -0.5 in to 0.5 m with the diaphragm at 
0.0 in is used in the calculations. Grids of 200, 400, 800, 1600 and 3200 cells are examined. All 
solutions arc computed at time equals 4.5 x 10 -4 sec. 

Results 

Results shown from the parametric study systematically progress toward obtaining the best 
time-accurate LAURA solution by isolating the individual parameters. The progression of the 
results examines the contributions of relaxation sweeps, Courant number, inviscid relaxation fac- 
tor, and grid spacing toward a time-accurate LAURA solution. The investigation begins with a 
comparison between first and second order accurate spatial differencing. 

Figure 1 displays the nondimensionalized density distribution, as calculated using first and 
second order accurate spatial differencing, compared to the exact solution. The Courant number 
for these computational solutions is 0.1. Twenty relaxation sweeps are taken for each time step over 
a computational grid of 200 cells. The second order accurate solution is computed in approximated 
2600 seconds on a Sun SPARCstation. Figures 2 and 3 show the pressure and velocity distributions 
for the comparison in Fig. 1. The additional accuracy obtained from second order differencing 
legitimizes its sole use in the remaining results of this study. Trends appearing in the velocity 
and pressure distributions are evident in the density distribution, notably, the dissipation existing 
across the shock and the expansion head. Thus, the results of the parametric study are discussed 
with regard to the density distributions of the solutions. 

The effects of the number of relaxation sweeps can be isolated by holding the Courant number 
constant. Figure 4 displays the solutions for a constant Courant number of 2 and relaxation sweeps 
of 10, 20, and 200. In this plot the number of relaxation sweeps is seen to affect the time-accuracy 
of a solution. Figure 5 examines the relaxation sweeps effect for a Courant number of 1. The only 
noticeable difference between 10 and 200 relaxation sweeps occurs across the shock; the 10 and 20 
relaxation sweeps solutions overlap. Figure 6 shows that for cn=0.5, 10 relaxation sweeps is enough 
to obtain a time-accurate solution. 

The two components of overall solution accuracy are time-accuracy and quality. Time-accuracy 
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is revealed by the calculated location of the shock, contact surface, and expansion wave compared 
to the locations predicted by the analytic solution. The quality is the degree of dissipation of 
the solution. Examining Figs. 4, 5, and 6 reveals that the number of relaxation sweeps affects 
the time-accuracy of a solution and the Courant number affects the dissipation of a solution. In 
Fig. 7 it is apparent that 5 relaxation sweeps is enough to obtain a converged solution for cn=0.1. 
Comparing the dissipation in Fig. 7 with that of Figs. 4, 5, and 6, it is clear that decreasing the 
Courant number produces solutions of less dissipation. Decreasing the Courant number any more 
than 0.1, however, increases the computational cost for a negligible gain in accuracy. 

If time-accuracy can be achieved for less than 5 relaxation sweeps by varying the inviscid 
relaxation factor, computational cost will be cut. Figures 8 and 9 display the effects of r/jjvv on 
the solution accuracy. In Fig. 8 the rfj^y is lowered from 2.0 to 1.5 for the 2 relaxation sweeps 
case. It appears that lowering the rfisv to 1.5 helps converge a solution that is not time-accurate. 
Figure 9 shows instability at the expansion tail caused by rf/NV= 1.5 for a Courant number of 2. 
This instability might be removed by updating the Jacobian of Eq. 7 more often than once per time 
step. The study also revealed, however, that a solution with r//wv=2.0 may actually converge to 
time-accuracy before a solution with rfjNv= 1-5 regardless of Courant number. 

The results thus far have used a computational grid of 200 cells. The highest quality LAURA 
solution is the case with cn=0.1, j //^v'= 2.0, and 5 relaxation sweeps. The remainder of the study 
will examine the effect of increasing the grid resolution. 

Figure 10 compares the 200, 400, and 800 grid cells solutions. It is apparent that increasing 
the grid resolution improves the quality of the solution. The price of this increased quality is 
computational cost. The solution for 800 cells is computed in 41000 sec. Figure 11 reveals that 
increasing the number of grid cells to 3200 produces a better approximation of the exact solution. 
At this scale there is little noticeable difference between the 800 grid cells solution and the 3200 
grid cells solution. Figure 12 displays in detail the shock and contact surface shown in Fig. 11. 

Increasing the number of grid cells significantly adds to the computational cost of the solutions. 
At increased grid resolution, it would be advantageous if an accurate solution could be found for less 
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than 5 relaxation sweeps. Figure 13 shows that across the shock a noticeable difference exists for the 
2 and 5 relaxation sweeps solutions. This result points out that while increasing the grid resolution 
significantly decreases the dissipation there is a less pronounced effect on the time-accuracy of a 
LAURA solution. 

Conclusions 


This parametric study of a time-accurate version of LAURA applied to inviscid, perfect-gas 
shock tube flow is an investigation into the trade-offs between overall solution accuracy and com- 
putational cost. The results of the study lead to some interesting conclusions concerning the effect 
of the parameters, the best LAURA solution, and computational cost of the solutions. 

The four parameters examined are Courant number, relaxation sweeps, inviscid relaxation 
factor, and grid spacing. The parametric study indicates that these coupled parameters have, in 
fact, isolated effects on the solutions. The Courant number and grid spacing significantly affect 
the dissipation of the solution, while the number of relaxation sweeps and the inviscid relaxation 
factor afTect the time-accuracy of the solution. As would be expected, second order accurate spatial 
differencing results in higher resolution quality than first order accurate spatial differencing. 

Finding the best time-accurate LAURA solution was a primary goal of this study. The best 
solution is defined as the solution with the highest degree of accuracy for the least computational 
cost. For a grid of 200 cells, the best solution is cn=0.1 with 5 relaxation sweeps. The results of 
the study indicate that LAURA is capable of producing extremely accurate solutions by increasing 
the number of grid cells. This study was limited to the maximum of 3200 grid cells in consideration 
of excessive computational time. It appears that by increasing the number of grid cells the quality 
of the solution will continue to improve. 

The price of highly resolved solutions is considerable computational cost. The solution of the 
cn=0.1, 5 relaxation sweeps, and 200 grid cells case is computed in about 43 minutes. Increasing 
the grid to 800 cells, while holding Courant number and number of relaxation sweeps constant, 
increases the computational time by a factor of 16. 
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Figure 1 . Effect of spatial-differencing accuracy on the nondimensional 

density distribution. 
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Figure 2. Effect of spatial-differencing accuracy on the 

nondimensional pressure distribution. 
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Figure 3. Effect of spatial-differencing accuracy on the 

nondimensional velocity distribution. 
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Figure 4. Effect of relaxation sweeps on the nondimensional 

density distribution for cn=2. 



o 



* 

CL 


13 


x,m 

Figure 5. Effect of relaxation sweeps on the nondimensional 

density distribution for cn=1 . 
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Figure 6. Effect of relaxation sweeps on the nondimensional 

density distribution for cn=0.5. 
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Figure 7. Effect of relaxation sweeps on the nondimensional 

density distribution for cn=0.1 . 
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Figure 8. Effect of the inviscid relaxation factor on the 
nondimensional density distribution for cn=0.1. 
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Figure 9. Effect of the inviscid relaxation factor on the 

nondimensional density distribution for cn=2 
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Figure 10. Effect of grid cell spacing on the nondimensional 

density distribution. 
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Figure 1 1 Effect of grid cell spacing on the nondimensional 

density distribution. 
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Figure 12. Effect of grid cell spacing on the nondimensional 

density distribution across the shock and contact surface. 
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Figure 13. Effect of relaxation sweeps on the 

nondimensional density distribution. 
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